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Abstract 

We apply a new formalism to derive the higher-order quantum kinetic expansion (QKE) for 
studying dissipative dynamics in a general quantum network coupled with an arbitrary thermal 
bath. The dynamics of system population is described by a time-convoluted kinetic equation, where 
the time-nonlocal rate kernel is systematically expanded on the order of off-diagonal elements 
of the system Hamiltonian. In the second order, the rate kernel recovers the expression of the 
noninteracting-blip approximation (NIBA) method. The higher-order corrections in the rate kernel 
account for the effects of the multi-site quantum coherence and the bath relaxation. In a quantum 
harmonic bath, the rate kernels of different orders are analytically derived. As demonstrated by 
four examples, the higher-order QKE can reliably predict quantum dissipative dynamics, comparing 
well with the hierarchic equation approach. More importantly, the higher-order rate kernels can 
distinguish and quantify distinct nontrivial quantum coherent effects, such as long-range energy 
transfer from quantum tunneling and quantum interference arising from the phase accumulation 
of interactions. 
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I. INTRODUCTION 



Quantum dissipation plays a key role in understanding quantum dynamic processes. The 
interaction between a quantum system and its surrounding environment causes an irre- 
versible loss of the energy and coherence of the quantum system. The relaxation and de- 
coherence times are the limiting factor of the quantum computation and quantum informa- 



tion 



lj . In the Caldeira-Leggett model, the change of the dissipation strength can interpret 



quantum tunneling and localization in macroscopic systems [2|, 



a. 



For many years, the 



solvent modulation in chemical reactions and quantum transport processes have attracted 
a lot of attentions (4, sj]. Incorporated with the description of the solvent reorganization, 
the Marcus theory is able to explain essential features of electron transfer 6| . In the recent 
two-dimensional (2D) electronic light spectroscopy, long-lived quantum coherence and wave- 
like dynamics have been found in natural light-harvesting protein complexes |3] and organic 
conjugated polymers [8|, which also triggers studies on the energy transfer optimization from 
the dissipation induced by the protein environment js-15|. To understand the nontrivial ef- 
fect of quantum coherence in the energy transfer, we need to study the underlying quantum 
dissipative dynamics beyond the conventional Forster resonance energy transfer (FRET) 
theory [16]. 

An accurate and reliable approach to compute quantum dissipative dynamics is a long- 
lasting but difficult theoretical problem. A large number of methods have been developed 
under many different frameworks, such as the Nakajima-Zwanzig projection operator txs| , 



the Feynman- Vernon influence functional approach [19j, and quantum stochastic noises for 



mulation 



20|422]. 



Redfield equation 



he second-order perturbation methods, such as Fermi's golden rule rate, 
23], generalized Bloch- Redfield equation js), [3], and noninteracting-blip 



approximation (NIBA) [3[, are derived in the limit of either a weak or strong system-bath 



interaction. In the variational polaron approach 



25] , a self-consistent reference can partially 



improve the prediction of the second-order perturbation. With a classical bath, the Ha 



ken- 



Strobl-Reineker (HSR) model 26N29j| and other quantum-classical mixed methods 30 



describe dissipative dynamics at high temperatures. The dissipative dynamics under a quan- 
tum harmonic bath can be evaluated by many sophisticated methods, such as the semiclassi- 
cal initial value representation (SC-IVR) 34J, the iterative linearized density matrix (ILDM) 



propagation 



35|, the quasi-adiabatic propagator path integral (QUAPI) 36|, the path in- 
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tegral Monte Carlo |37|, |38], etc. Despite their successes, these methods can be numerically 
expensive and becomes difficult for long-time dynamics. If the time correlation of the har- 
monic bath can be expanded as a sum of exponentially decaying functions, the hierarchy 
equation approach can accurately predict quantum dissipative dynamics by expanding over 



auxiliary fields 39H43] . However, the hierarchy equation is numerically difficult for large- 



scale systems and the bath with a long-tail correlation, and converges slowly for strong 
system-bath interactions and low temperatures. 

Hopping kinetics of the Fermi's golden rule rate is often considered as a 'classical' descrip- 
tion of quantum dissipative dynamics, although the two-'site' quantum coherence is included 
in the rate expression [lllQj. The interesting quantum phenomena beyond the second-order 
hopping kinetics can be attributed to nontrivial quantum effects of multi-'site' coherence, 
e.g., long-range transfer (tunneling) and quantum interference. The temporal correlation of 
bath is also crucial in understanding the full quantum dynamics. The higher-order bath re- 
laxation effect is caused by the deviation from the system-bath factorized reference state. In 
addition, the second-order hopping kinetics cannot predict the detailed balance of quantum 
dynamics, i.e., the Boltzmann equilibrium distribution including both the system and the 



bath 



44 



45]. The comparison of the second-order hopping kinetics and the full quantum 
dynamics in our previous paper jll| has revealed the integrated behavior of the above effects, 
together with the flux network analysis. To distinguish and quantify nontrivial quantum 
effects, we need a systematic expansion procedure to calculate every higher-order correction 
beyond the second-order hopping kinetics, which is almost impossible in many sophisticated 
theoretical methods discussed above. Following the stationary approximation for the co- 
herent term, the kinetic mapping of quantum dynamics allows us to identify the multi-site 



29]. However, the theoretical method 



quantum coherence term by term in the HSR model 
for a general quantum bath is still missing. 

To address the above concerns, we will apply a general non-Markovian quantum kinetic 
equation, where the time-nonlocal rate kernel is obtained by a systematic expansion ap- 
proach. This higher-order quantum kinetic expansion (QKE) method presents a rigorous 
mapping from a quantum network to a kinetic network, which helps us to identify nontrivial 
quantum effects and bath relaxation beyond the traditional classical description. In addi- 
tion, the higher-order QKE is expected to serve as a numerically reliable method to calculate 
the population dynamics. For simplicity, we focus on a product state between the system 
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and the bath at zero time, and assume that quantum system is initially prepared in the 
population subspace without quantum coherence. This initial incoherent assumption is ac- 
ceptable for an initially localized (quasi-) particle in the quantum transport process, e.g., an 
exicton after absorbing incoherent sunlight in a natural light-harvesting protein complex (4]. 
Accordingly, the bath induced fluctuation is assumed on the diagonal elements of the system 
Hamiltonian. Since the system is not defined in its eigenbasis (the system Hamiltonian is 
not diagonal), the diagonal fluctuation can still lead to both relaxation and decoherence. As 
we discussed, many theoretical approaches require the presumption of the harmonic bath, 
which is often considered as a good approximation under many circumstances. In complex 
environments such as a protein backbone, anharmonicity however could be relevant even in 
a fast energy transfer process. Here, the higher-order QKE is free of the harmonic bath pre- 
sumption, although the additional numerical implementation is required in the anharmonic 
bath. 

An essential part of our theory is the systematic expansion of time- nonlocal kinetic rate 
kernels. For the two-site system, the bath relaxation effect was calculated in electron transfer 
following a term-by-term comparison in integrated population between microscopic expan- 



sion of quantum dynamics and a formally exact non-Markovian kinetic equation 



46|. Here 



we will extend this procedures to a multi-site system. Although some of the higher-order cor- 



rections have been studied in the past 



47H50|. the higher-order QKE in this paper provides 



a systematic formalism of obtaining the general expression of the time-nonlocal rate kernels 
and unify the bath relaxation and the multi-site coherence under the same framework for a 
multi-site quantum network. 

The paper is organized as follows: In Sec. Ull we will integrate the time evolution of 
quantum coherence and project the Liouville equation to be a closed dynamic equation of 
system population. In Sees. IIIH we will develop the non-Markovian higher-order quantum 
kinetic expansion method. The time-nonlocal kinetic rate kernel is generalized from the 
simplest second-order, i.e., the NIBA expression, to an arbitrary fc-th order. The time 
correlation function formalism in the Hilbert space provides rigorous expressions of rate 
kernels in an arbitrary bath. In Sec. IIVt we will focus on the harmonic bath and apply 
the displacement operator and the cumulant expansion to derive the analytical expressions 
of rate kernels. In Sec. |V] the higher-order QKE is applied to four model systems for its 
reliability. In addition to its numerical accuracy, we identify the bath-induced slow-down in 
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the quantum transport rate, and two nontrivial quantum coherent effects, quantum tunneling 
and quantum phase interference. In Sec. IVIt we will conclude and discuss the higher-order 
QKE method. 



II. POPULATION DYNAMICS PROJECTED FROM THE LIOUVILLE EQUA- 
TION 

For an arbitrary open quantum network, the total Hamiltonian is written as H = 
Hs + Hb + Hsbi where H$ and Hb denote the bare system Hamiltonian and the bath 
Hamiltonian, respectively. The interaction between the system and the bath is described 
by Hsb- The bare system Hamiltonian Hs is defined in its ^-dimensional (iV-D) Hilbert 
space. In the single-excitation manifold of a Frenkel exiciton system, the n-th basis of the 
Hilbert space, \n) = |0, • • ■ , 1, 0, ■ • • ), represents a combination of one excitation state at 
the n-th local chromophore site and the ground state at all the other chromophore sites {4]. 
The total Hamiltonian H can be also expanded in the iV-D system Hilbert space using 
\n)\b) = |0, bi, ■ ■ ■ ; 1, b n \ 0, b n+1 ; ■■■), where \b n ) = \b nl , b n2 , - - • , b nMn ) is the complete basis 
set of M n {— > 00) bath modes sorrounding |n). Here we consider a bath-induced fluctuation, 
HsB;n = J2bb> HsB;n,b;n,b'\b) (b'\, over each diagonal element (e n ) of Hs- The fluctuation over 
off-diagonal elements ( J mn ) of Hs is not included in our current model; this approximation 
is often applied in the study of energy and charge transfer . The total Hamiltonian is 
given as 



where H n = e n + Hb + Hsb-,u is implicitly a quantum operator of bath. 

The time evolution of the total density matrix p(t) for the system-bath Hamiltonian is 
governed by the Liouville equation, p = —iCp, where £ = [H, • • • ] is the Liouville super- 
operator. The Planck constant h is treated as a unit throughout this paper. With respect 
to the system basis {|^)}, we divide p into two sets: population p P = {p nn } and coherence 
Pc = {Pmn{^m)}- Each element, p nn or p mn , is a quantum operator of bath and implicitly 
includes information of the entangled system and bath, i.e., p nn = J2bb> Pn,b;n,v\b){b'\ and 
Pmn = Ylibb' Pm,b;n,b'\b)(b'\. In this paper, we will derive our theory using both Hilbert and 
Liouville frameworks. To distinguish notations in these two frameworks, we will use 'state' 




(1) 
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to specify a density state (population and coherence) in the Liouville space, unless otherwise 
explained. The wavefunction basis of the Hilbert space will be referred as 'site', consistent 
with the single-excitation manifold in the multi-site exciton network. The orginal Liouville 
equation is divided into two coupled equations, 

p F = -i£ F p F - iCpcPc, (2a) 
Pc = ~i£cpc ~ ^CcpPp, ( 2b ) 

where the two subscripts P and C denote population and coherence of the system, respec- 
tively. The total Liouville superoperator is expressed in a block matrix form, 



C= ~* ~~ . (3) 




Next the coherence vector is integrated to yield 

p c (t) =Mc(*)pc(0) -i I dr 2 ^r 1 5 ri+r2i fW c (r 2 )£cppp(n), (4) 

Jo Jo 

where Uc{t) = exp(-iCct) is the time evolution matrix of coherence in the Liouville space. 
For simplicity, we assume zero initial coherence, pc(0) = 0, so that the first term on the 
right hand side of Eq. fl4]) vanishes. More complicated initial conditions will be left in 
the future. In the two-'site' system, £ c and Uc are diagonal in the system basis, i.e., 
£c ; i2,2i = and Wc ; i2,2i = 0. In the iV(> 2)-'site' systems, £c is no longer diagonal, but 
diagonal and off-diagonal elements might be distinguished by their orders of magnitude, e.g., 
\£>C;mn,mn\ 3> \£c;mn,m'ri{^mn)\, m the strong damping limit. We express £c as a sum of the 
diagonal matrix C§. mn m , n , = C c - mn ^ mn 5 m ^ m 5 n > >n , and the remaining term, C§ = C c - C£ \ 
As shown in Sec. IIII D\ the separation of and Cq is equivalent to the separation of H n 
and J mn . 

Expanding the coherence vector pc in the order of and substituting the result into 
Eq. (j2ap . we obtain a closed time evolution equation of population, 

oo „ 

fr(t) = -iCrpAt) + J2 (-i) k £pcU { c\T k )££ ) U { °\T k - 1 ) ■ ■ ■ C^U^ir^CcPPpin). (5) 

k=2 J 

Equation (jSJ) shows that a population flow always passes intermediate quantum coherence 
states in the Liouville space of density states. In the two-site system, population and co- 
herence states appear subsequently, since the direct interconversion is not allowed between 
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the two coherence states, pn and p 2 i (£c ; i2,2i = 0) [3j. In the N(> 2)-site system, the 
direct interconversion between different coherence states is allowed by a nonzero 'interac- 
tion', Cfj, from the multi-site quantum coherence. The interactions responsible for the 
transition between coherence and population states are £pc and £cp- All the three terms, 
Cpci £cp), ar i se from the off-diagonal elements J mn of the bare system Hamiltonian, 
and will be counted together in the expansion order of the final quantum kinetic equation 
(QKE). For conciseness, we introduce the fc-th order time-nonlocal population transition 
matrix, 

W« = -(-^^^(^^^(rfc-x)---^^^)^, (6) 

where k is the total number of £ terms, including (£q , Cpc, £cp)- Using the complete 
transition matrix, W = + + W^ 4 * 1 + ■ • • , we formally rewrite the time evolution 
equation of population as 

pp{t) = -zXppp(t) - W*p P , (7) 

where the symbol * represents a general time convolution form, 

X *Y = dr x - ■ ■ dTidT i+1 ■ ■ ■ drjXfa, ■■■ , Ti)Y(r i+1 , ■■■ , T j )5 n+ ... +n+Ti+1+ ... +Tjit , (8) 
Jo 

for two arbitrary functions X(t\, ■ ■ ■ ,T{) and Yfc+i, • ■ ■ , Tj) of time. Equation (JTj) is equiva- 
lent to the projection of the original Liouville equation onto the population subspace without 
averaging over bath. A pratical computation of the reduced system population dynamics 
relies on further simplifications introduced in Sec. IIHI 



III. NON-MARKOVIAN HIGHER-ORDER QUANTUM KINETIC EQUATION 

A. Local Born Approximation and Multi-Site Quantum Coherence 

In microscopic quantum systems, a useful physical observation is the time scale sepa- 
ration between different degrees of freedom. In this subsection, we assume that each lo- 
cal bath b n can instantaneously relax to its Boltzmann equilibrium density state, = 
exp(— /3if n )/Trb{exp(— /3H n )}, at any moment. Thus, the n-th element of the transient to- 
tal population vector is written in a product form, pp ]n (t) = P n (t)p e b q 



46|. This locally 
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fast bath (Born) approximation is different from p(t) = ps{t)pl q and p£ q oc exp(—PHs) 



in the Redfield equation 23]. The time scale separation in many realistic systems is not 
always satisfied, so that we will relax the above local Born approximation in Sec. IIIIBI and 
systematically include the contribution of bath relaxation. 

In this subsection, we discuss higher-order corrections of system coherence, i.e., multi- 



site coherence 



11 



29[, under the local Born approximation. After averaging over bath in 



Eq. ([7]), the time evolution equation of the system population is given by, 

P(t) = -Tr 6 { [W (2) + W (3) + W (4) + •••]* p c B q } * P, (9) 

where p 1 ^ 1 is the vector of the equilibrium bath density state at each local site, i.e., p 1 ^ 1 = 
(Pfei 5 Pfe' ' ' ' ) T- Here the superscript T denotes matrix transpose, and the symbol * defines a 
matrix product, (X -kY) mn = X mn Y n , of a matrix X and a vector Y. To be consistent, each 
local bath is in equilibrium initially, i.e., pp ;n (0) = P n (0)p^ q , which belongs to the class-B 



preparation in Ref. 51|. 



For conciseness, we introduce two quantum bath operators, (X = Tr^jX, and X) = 
X * p C g}- Consequently, we define the bath average, (X) = Trj, {X * p^}, and the projection 
onto the bath equilibrium density state, ) ( = *p^} Tr^ { . Equation fl9]) is simplified to 
a time-convolution form, P(t) = — (W) * P. Compared to the second-order truncation 
approaches such as Fermi's golden rule rate and NIBA, Eq. (J2J) systematically includes 
time-nonlocal corrections of multi-site quantum coherence, {W^, W^ 4 **, ■ ■ • }, resulted 
from direct interconversion of coherence in the Liouville space. 

B. Bath Relaxation Effect 

Normally, the bath requires a characteristic time scale to adjust to the system change 
and relax back to equilibrium. The resulting memory kernel can be crucial for long-lived 



quantum coherence in light-harvesting systems 



reactions Therefore, we need a systematic and reliable way to include the contribution 
of bath relaxation beyond the Born approximation. To compute higher-order corrections 
from both quantum coherence and bath relaxation, we extend an approach previously for 



and solvent-modified electron transfer 



the two-site system 46J to the general TV-site system. 



Integrating the time differential equation in Eq. (J7|), the total population vector in a non- 
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equilibrium bath is written explicitly as a time-convolution form, pp(t) = Wp(£)pp(0) — Up * 
W * pp, where ZYp(t) is the time evolution matrix of population, Up(t) = exp(— iCpt). For 
a product initial state, pp ;n (0) = P n (0)p|j q , the bath average leads to the system population 
in the form of 

P{t) = P(0) - [*(W (2) >*] P(0) - [*(W (3) >*] P{0) 

+ [* ((W (2) * Up * W (2) ) - (W (4) )) *] P(0) + • • • , (10) 

where the expansion order is the total number of the Liouville superoperators, including 
Cq, Cpc, and Ccp- Similar to that in Ref. 46], the notation of [*X*] defines a time 



convolution with the unit function in both the first (n) and final (r^) time steps, i.e., 

[*X*] = / dridr 2 • • ■ dT k X(r 2 , ■■■ , T k „ 1 )5 Tl+T2+ ... +Tktt . (11) 
Jo 

On the other hand, we can formally assign a time-nonlocal kinetic equation, 

p(t) = _/C * P, (12) 

to describe the time evolution of system population P, where /C is the time-nonlocal quantum 
rate kernel. Similar to W, the rate kernel fC can be expanded as K = /C^ + + /C( 4 ) + • • ■ , 
in the order of the £ terms. The integration of Eq. (TT2|) then leads to 



P(t) = P(0) - [*/C (2) *] P(0) - [*/C (3) *] P(0) 

+ [* (K® * - K^) *] P(0) + ■ ■ • . (13) 

The term-by-term comparison between Eqs. (fTOl) and (fT3|) determines the explicit forms of 
the quantum rate kernels, e.g., 

£(2) = ^(2)^ ( 14a) 

^(3) = ^(3)^ (Mb) 

/CW = (W (4) ) - [(W (2) (T- 4 )Wp(r 3 )W (2) (r 2 )) -/C^(r 4 )/C^(r 2 )] . (14c) 
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C. Higher-Order Quantum Rate Kernels and Kinetic Mapping 



Extending the procedure in the previous subsection to higher orders, we can straightfor- 
wardly derive the general form of the fc-th rate kernel, given by 

/C (fc>3) (r 2 , • • • , r k ) = (#)) - [{W [kl) UpW^) - tf*Otf*»>] 

fcl,fc 2 >2 

+ 6 k 1+ k 2+ k 3 ,k [{W {kl) UpW {k2) U P W [ki) ) - /C( fel )/C (fe2) /C (fe) ] + ■ • • . (15) 

The right hand side of Eq. (TIB]) is terminated when each index ki of /C^^/C^ 2 ^ • • - KS k ^ ■ ■ ■ 
in the final summation term is equal to either 2 or 3. The summation terms subsequently 
changes between positive and negative signs. The time variable sequence follows the same or- 
dering, {r 2 , • • • , Tfe_i, T/c}, in each summation term. Comparing Eq. (fT5|) with the expression 
in Sec. IHI Al we observe that in addition to the correction (W^) from multi-site quantum 
coherence, the bath relaxation (system-bath entanglement) also influences quantum dynam- 
ics due to the fluctuation around the reference density state of the local Born approximation 
(the local equilibrium density state of bath). For example, the second term of K,^ can be 
simplified to (W^SUpW^), where Slip = Up—){ represents the deviation from the local 



43], Eq. ([15]) includes the odd-fc 



Born approximation. Compared to the expression in Ref. 
terms from the imaginary accumulated phases. 

The non-Markovian quantum kinetic equation in Eq. ( Tl2l) together with kinetic rate 
kernels in Eqs. fll4l) - fll5l) constructs a rigorous theoretical framework, i.e., the higher-order 
quantum kinetic expansion (QKE), to compute the time evolution of the system population 
in a quantum network. The quantum dynamics of the density matrix in the N 2 -D Liouville 
space is mapped onto kinetics in the iV-D population space. In the HSR model, such kinetic 



mapping was developed based on the stationary approximation of coherence 29|- In our 



current formalism, kinetic mapping is generalized for the arbitrary N-D quantum network 



using the non-Makrovian rate kerne 
rate kernel in the NIBA approach 



fC. The leading-order expansion, /C( 2 \isthe same as the 
2^. The time integration of K,^ recovers Fermi's golden 
rule rate, which is often considered as the 'classical' description of kinetics. As corrections 
to K,^ 2 \ higher-order rate kernels Kt® can identify and quantify various nontrivial quantum 
coherent effects, which will be demonstrated by examples in Sec. [V] 
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D. Quantum Kinetic Rate Kernels Expressed in Hilbert Space 



To compute quantum kinetic rate kernels, we express the superoperators C and U{t) as 
functions of the Hamiltonian H and the time evolution operator U(t) in the Hilbert space. 
Based on its definition, CX = [H,X], the Liouville superoperator £ is expanded to be, 
C mnM X k i = H mk X k i5 n j - X kl Hi n 5 m , k , where the positions of H mk and H in are usually fixed 
since these two Hamiltonian elements can be quantum operators of bath, the same for X. 
In this paper, we ignore the fluctuation around off-diagonal Hamiltonian elements, which 
leads to the following scalar forms, 

£pC;m,M = Jrnk&rn,l ~ Jlm&m,ki (16a) 
£cP;W,m = Jkmh.m ~ Jml$k,m, (16b) 

^C-MhMh = Jk^h-h - Jhhhuki- ( 16c ) 

The other two Liouville superoperators, C-p and Cq , are diagonal in the system basis. Each 
diagonal element of the two corresponding time evolution matrices behaves as 

Up;m(t)X = U m (t)XU+(t), (17a) 
4° ; L(*)* = U n (t)XU+(t), (17b) 

where U n (t) = exp(—iH n t) is the time evolution operator of the local site basis \n, b) in the 
Hilbert space. 

Next we substitute Eqs. (TT6]) and ffTTl) into the expressions of quantum kinetic rate kernels 
derived in Sees. IIIIBI and IIII CI For example, the second-order kinetic rate kernel becomes 

K%U*n) = -2|^n| 2 Re Tr 6 {U+Wn(t)&} , (18) 

where 'Re' denotes the real part of a complex variable and the imaginary symbol 'Im' will 
also be used in this paper. The higher-order kinetic rate kernels can be similarly obtained. 
So far our derivation is rigorous and general: The surrounding bath can be an ensemble 
of harmonic or anharmonic oscillators with an arbitrary spectral density. The bath can 
alternatively be defined by nuclear motion of molecules and atoms, following quantum or 
a classical dynamics. The system-bath coupling Hsb can follow any functional forms in 
addition to the regular bilinear form. For complex baths, numerical simulation will be 
required to calculate the rate kernel of different orders. 
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IV. HIGHER-ORDER QUANTUM KINETIC EXPANSION FOR A HARMONIC 
BATH 



In the remainder of this paper, we will focus on a blinear coupling between the system and 
a harmonic (Boson) bath. In this section, we will derive analytical expressions of rate kernels 
required in the higher-order QKE. With the creation (af) and annihilation (oj) operators 
for the i-th harmonic oscillator, the diagonal Hamiltonian element for each system site \n) 
is written as 

H n = E n + 22 U i a t a i + UiXni ( ai + ( 19 ) 

i i 

where the quantum zero-point energy Ui/2 is ignored and the coefficient x n % is the system- 
bath coupling strength reduced by the frequency of the i-th harmonic oscillator. Quantum 
operators of different harmonic oscillators are assumed to commute with each other, i.e., 
[a^ + \ <4 ] = for i ^ j. Equation (ITTJ]) implicitly assumes a universal environment for all 
the system sites, and an alternative approach is to apply an isolated environment for each 
site; these two methods can lead to the same result. 



A. Canonical Transformation of the Displacement Operator 

The trace of time-dependent operators over the quantum harmonic bath can be solved 



by many theoretical techniques, e.g., the path- integral method [19l. |51|. Here we will apply 
a canonical transformation method together with the cumulant expansion. 

The displacement operator, G n = exp \Y^ i x ni (a~l — a*)] , is used to diagonalize the bath- 
modulated diagonal Hamiltonian element, resulting in system-bath decomposition, 

GnHnG" 1 = e n + H B , (20) 

where a shift appears in the diagonal energy, e n = e n — Ylii^^ni- Although the same 
canonical operator is applied in the polaron method, our general quantum kinetic equation 
formalism does not rely on the concept of polaron, as stated in Sec. IIHI The diagonalization 
in Eq. (120]) allows us to factorize the local time evolution operator into a product form, 

U n {t) =U S ,n{t)[G; i l U h {t)Gnl (21) 
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and express the local bath equilibrium state operator as 

pZ = G~ l pTG n) (22) 

where Us- n {t) = exp(—ie n t) is the time evolution function of the displaced system, and the 
other two operators, U b {t) = exp(— iif^t) and p e b q oc exp(— (3Hb), only depend on bath. In 
the Heisenberg picture, the time-dependent displacement operator is then written as 

Gn{t) = U+{t)G n U b {t) = exp |^x m [a+(t) - <n(t)]\ , (23) 

where ai(t) = a^e -1 ^* and af (t) = a^e luJit are time-dependent annihilation and creation 
operators, respectively. Substituting the above results into the second-order quantum kinetic 
rate kernel K^ 2 \ we arrive at 

JC { l ¥m) (t) = -2|J mn | 2 Re e^Tr {G nm (t)G mn p?} 

= 2|</ mri | Re e (G nm (t)G mri ) b , (24) 

with G mn G m G n exp \ x mrit i(b^ ^«)] > •^mn,i %mi x n i, and s nm E n E m . The 
average, (X) b = Tr b {Xp e b q }, is taken over the decoupled bath. By extending this method 
to higher-order expressions, we observe that all the quantum kinetic rate kernels are fully 
determined by multi-time correlation functions of the canonical operator G mn . 



B. Time Correlation Functions of Position Shift Operator 

In this subsection, we derive the general form of multi-time correlation functions of G mn 
in the harmonic bath. 

Applying the quantum thermal average of the harmonic bath, we obtain the analytical 
form of the two-time correlation function, 

(G mini (t)G mini (0)) b = exp{ 

5^712)12, mini 

(*)}, (25) 

where 

SWa.mmi (t) = ^2 Xrn^^iXmim.i [(1 - COS CJjt) COth(/3Wi/2) + I sin Uit] . (26) 
i 

In practice, we can assume a 'spatial' correlation, x m iX n i = c mn x 2 , between each pair of 
system sites, \m) and In the spin-boson model, a perfect negative correlation, c mn = 

13 



2£mn — 1, can be deduced from the Pauli matrix a 



jsjl. In energy transfer systems, a 
zero spatial correlation, c mn = 5 m , n , is often used Qj. For a continuous bath, the spectral 
density is defined by J(u) = ]\ Ji u 2 xf5(u — Ui), and Eq. (|2T)j) is simplified to g m2 n. 2 ,m 1 ,n 1 if) = 

^m\ni,m2ni9\^) with 



g(t) = / gL>[J(w)/(X> 2 ] [(1 — coswt) coth(/3o;/2) + i sin cut] . 
Jo 



(27a) 
(27b) 



The above procedure can be straightforwardly extended to multi-time correlation functions 
following the cumulant expansion of the Gaussian distributed noise. In general, the fcth-order 
time correlation function reads 



[Gm k n k {tk)Gm k _ 1 n k _ 1 (*fc-l) ' ' •G f mmi(*l)) 6 — exp 



. 3=2 j'=l 



,(28) 



where the index set, {m^m^-i, • • ■ , mi}, is the permutation of the original set of indices, 
{njunk-i, ■ ■ • ,ni}. An additional constraint, = m, is needed to close the index loop, as 
required by taking the trace. 



C. Three Leading-Order Quantum Rate Kernels in the Harmonic Bath 

Through a tedious but straightforward derivation, we obtain the analytical forms of 
quantum rate kernels in the harmonic limit. Here we summarize and discuss the result of 
the three leading order rate kernels, which will be applied to examples in Sec. [V] 

In the second-order quantum rate kernel, i.e., the NIBA rate kernel, each off-diagonal 
element is written as 

^mi(^m)( r 2) = -2| J mn | 2 Reexp {-[ie r . mn r 2 + s mn 5f(r 2 )]} , (29) 

with Sm n = s mn ,mn = 2(1 — c mn ). The diagonal element /Cnn is calculated by a summation, 
KXn = — Ylirn{^n) ^rnn- Following the original equation of the total density matrix in Eq. (JHJ), 
we can interpret each term of the time convolution, K,*P, as a dynamic trajectory of density 
states in the Liouville space, which determines the population evolution of system at the next 
moment. The diagrammatic representation of dynamic trajectories can clarify the effects 
of quantum coherence and bath relaxation in each term of /C. Figure [IK presents such a 
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dynamic transition, P n p^ n ~^ Pmn ~~ >* PmPWn accounted in Kmn- Here different circles of 
density states (population and coherence) are connected by arrowed lines, representing the 
direction of the dynamic transition. Each arrowed line is associated with a coupling J, 
which is the interaction responsible for the transition from one density state to the next one. 



Our diagrammatic representation resembles the pathways in Ref [47]] , but emphasizes the 
topology of the system Hamiltonian so that it is closer to kinetic mapping representation in 
Ref. [29| and easier to extract different dynamic behaviors in terms of expansion order. In 
addition, the factorized and unfactorized population states are plotted together to highlight 
the reduced kinetics in the population subspace. 

The higher-order rate kernels are corrections to K,^ 2 \ related to the multi-site quantum 
coherence and the bath relaxation. In detail, the third-order rate kernel is given by 



kSW^,t 3 ) = 2Im{j m „J nfc J fcm e j(£ - T2+ ^ T3) - F 3: fc ™ 



J nm Jmk Jkn 



e i{£ n kT2+e mk T 3 )-F 3 mkn _|_ e i{e nk T 2 +e nm T 3 )-F 3mnk 



(30) 



with F± abc = s catCb g(T 2 ) + s ab , ac g(±T 3 ) + s ba>bc g(r 2 + r 3 ). The summation over the extra 
system basis index k is implied in Eq. ( 130|) . and the same notation is applied to the other 
higher-order rate kernels. A typical dynamic transition, P n p b q n — > Pnm — > pkm — > PmP b q m , 
from the RHS of Eq. (1301) is plotted in Fig. [lb. The nonzero prefactor, JmnJnkJkmt requires 
a closed interaction loop in the system, so that /C® does not appear in a 1-D chain model 



under the nearest neighbor interaction 47|, [50]. For complex interactions, quantum phase 



interference can be significant in /C®, which will be demonstrated by an example of the 
three-site system in Sec. IV Dl 

The fourth-order quantum rate kernel can be divided into two terms depending on the 
time evolution operator in the intermediate step: /C^ th due to the bath relaxation (5Up) and 
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(i) 



due to the multi-site coherence (U^)- The ^ TS ^ t erm ^bath * s explicitly written as 



coherence 



(1) 



K (4) 

mn(^m);bath 



2Re ^ |t/ mn | It/mfcl 



S rim, km "4. A 1 



_|_ e *(£nmT2+<?mfc' r 4) g - [s mn g(r 2 )+s mfc g(r 4 )] ^s nmkm F 4A _ ^ 



I <^ mn I I Jnk I 



^mn,kn " 4A 1 



17 I 2 I T |2 



e i(£nkT2+£mkT-4) e -[s n kg(T2) + S mk g{-T4)] ^Q—S n k,mk F 4,A — \ \ 



with = flf(±r 3 ) - 5 , (r 2 + r 3 ) - #(±(t 3 + r 4 )) + g{r 2 + r 3 + r 4 ). As shown in Fig. El 
the dynamic transitions in /C^ th can be categorized into three types of diagrams: a) The 
interaction prefactor is | J mn | 4 , and the dynamic transition is within the sub-Liouville space 
of the starting and ending system sites, |n) and \m). A typical transitions is P n p^ n 
Pmn — > Pn — > Pmn PmPb^ m i where an intermediate population fluctuation occurs at site n 
because the non-equilibrium bath is entangled with the system, b) The interaction prefactor 
is \Jmn\ 2 \Jm{n)k\ 2 ■ In the dynamic transition, a coherent state between \m) (\n)) and an 
additional site \k) |m), |n}) is involved but the intermediate population fluctuation is 
still caused by the bath entangled with site \m) or \n). A typical transition is P n pl q n 
Pmn Pm — >• Pmk PmpVm- c ) The interaction prefactor is |«/ m fc| 2 |^nfc| 2 , so that the 
intermediate population fluctuation is caused by the bath entangled with the additional site 
| A;). A typical transition is P n p eq n — > Pnk — > Pj^~ ► Pkm — > PmPb^m- m ^ ne two-site system, 
only the first type of trajectories can appear 46|. In the A-site system, the bath relaxation 
can also induce a long-range transport from the second and third types of trajectories in 
Figs. |2b and c, in addition to the multi-site coherence. 

The other fourth-order term, Adherence' due ^° the m ulti-site coherence is explicitly given 

by 



(4) 

coherence;mn(7^m) 



2Re / 7 , 7, , T, T P i {^ m T2+S km T 3 +ei m T A )-F AB kmnl 

4- T uTui L T J{£mT2+e n kT3+£ nm T4)-F+ B m 
~ ° mk° Kl° ln° nm 

_A_ e i(£nlT2+E m lT3+E mk T4)-F- c . mlnk _|_ g i (e nl T 2 +E nk T 3 +E mk T 4 ) -F£ C ;knl™ 
— 7,7, 7,7, p i (£nkT2+E n mT 3 +Sl m T4)-F+ c fe; 

° mk° kn° nl° Im ° 

_|_ ^i(E nk T2+El k T 3 +El m T4)-F~ c . lknm _|_ R i{E nk T2+El k T 3 +£ mk T4)-F 4 ~ B . lkr 



,(32) 
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with 



F 4B;abcd = S ac,bc9{ T 2) + S acM g{±T Z ) + S adM g{±T A ) 

+s a d,cbg(r 2 + r 3 ) + s aCidb g(±(r 3 + r 4 )) + s bcM g(r 2 + r 3 + r 4 ), (33a) 

F 4C;abcd = S ac,bc9( T 2) + Sbd,ca9(±T 3 ) + S 6d;ad #(=FT 4 ) 

+Sbd,bcg(r 2 + r 3 ) + s dajCa 5f(±(r 3 + r 4 )) + s da)bc g(r 2 + r 3 + r 4 ). (33b) 

Notice that the site indices, m\ and 711, for an arbitrary oscillation frequency, £ mmi , cannot 
be identical in Eq. (|32|) . Based on the number of additional system sites in /C^ erence , we 
identify two types of multi-site coherence behaviors in the fourth order, and each type is 
further divided into two transition structures. As shown in Fig. [3^ and b, the first type 
°f ^coLrence-mn involves one additional system site: a) The site k interacts with either the 
starting (\n)) or the ending (\m)) site, and the interaction prefactor is \J m n\ 2 \Jn(m)k\ 2 - One 
example dynamic transition is P n p^ n -> Pmn -> Pmk -> p mn -> PmPl%,- b) The site 
interacts with both |m) and |n), and the interaction prefactor is | J m k\ 2 \Jnk\ 2 ■ One example 
transition is PnP b q n ~^ Pkn — ► Pmn Pmk PmP^ m - As shown in Fig- Efc an d d, the second 
type of A^cohercnce mn involves two additional system sites, \k) and \l), which interact with 
both \m) and \n) and form a closed loop: c) The starting and ending sites, \m) and \n), 
are interacted. One example transition is PnP b q n — > Pmn — > Pmk — > Pmi — > P-mP^m with the 
interaction prefactor J mn JnkJkiJim- d) The two sites, \m) and |n), are not interacted. One 
example transition is P n pt q n Pkn Pmn Pmi ~^ PmP c b q m with the interaction prefactor 
JmkJknJimJni- The long-range quantum transport in the linear chain system is explained by 
the first type of trajectories (4^, 0], whereas the four-site quantum interference is described 
by the second type of trajectories. 

V. EXAMPLES OF THE HIGHER-ORDER QUANTUM KINETIC EQUATION 

In this section, we apply the higher-order QKE to four model systems, examining its 
validity and reliability. To reduce the computation cost, we introduce the Markovian ap- 
proximation in the rate kernels. The time evolution of the population of system is changed 
to P — —KP, where K = + K^> + + • • • is the effective rate matrix defined by 
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the time integration of the rate kernel, 

poo J 1 

K^= / I^/C«(r 2 ,r 3 ,--- ,r k ). (34) 

^0 i=2 

The second-order effective rate, K^ 2 \ recovers Fermi's golden rule rate. The Markovian ap- 
proximation ignores the short-time quantum oscillation but can reliably describe the overall 
population dynamics. 



A. Kinetic Mapping in the Haken-Strobl-Reineker Model 



The first example is the Haken-Strobl-Reineker (HSR) model 26M29I] . where the bath is 
a Gaussian classical white noise. Without bath relaxation terms, only multi-site quantum 
coherence contributes to higher-order corrections, i.e., Kt® = (W^). A stationary coherence 



approximation was applied to derive kinetic mapping of the HSR model 29|. Here we will 
demonstrate that higher-order QKE leads to the exactly same result. 

The spectral density of white noise, J(uj) = r(3u/2ir, together with the high-temperature 
approximation, coth(/3o;/2) m 2/f3u, yields the time correlation function, g{t) = T\t\/2, 
where the imaginary part is omitted under the consideration of the classical noise. Without 
spatial correlation, c mn = 5 m ,m the second-order kinetic rate (i.e., Fermi's golden rule rate) 
is obtained as 

K (2) _ K (2) — 17 12 ZL ™n 

lv mn ll nm \ u mn\ -p 2 _i_ ?2 ' \°°) 



which is the same as that derived in Ref. 



29] 



All the higher-order corrections can be straightforwardly calculated by substituting the 
linear function of g(t) into expressions of K^ k \ To demonstrate the validity, we examine 
the closed-looped three-site model. Following Eqs. ( 130]) and fl34|) . the third-order correction 
from site 2 to site 1 is given by 

K® = -2 Im( Jl , 3J l 2J21 + J ! 3J l 2 ' 721 + Jl3J32J21 l , (36) 

I r2ir3i r32T3i r32ri2 J 



where T mn = T mn + ie mp . is the complex dephasing rate. Equation ( 1361) is identical to the 



result derived in Ref. 



291 ]. and the same conclusion is applied to all the other HSR systems. 
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B. The Bath Relaxation Effect in a Two-Site System 



The second example is a two-site system in a quantum bath (see Fig. Hk). This model is 
widely applied in the study of quantum transport and quantum phase transition. Without 
the coherence- coherence transition {Cq = 0), all the terms with VV ( - fc>2 ^ disappear and 
higher-order corrections only arise from the bath relaxation effect, differing from the pure 



46j , The bath relaxation effect in the 



multi-site coherence effect in the HSR model. In Ref. 
spin-boson model has been studied following the short-time asymptotic expression of g(t). 
We will extend the calculation to the donor-acceptor pair with the zero spatial correlation, 
Cmn — $mn- To compare with the exact quantum dynamics, we consider a quantum bath 
described by the Debye spectral density, which can be alternatively solved by the hierarchy 



equation 



33-142 1 . 



(37) 



The Debye spectral density is given by 

J(uj) = — —„ T 

where Q(u) is the Heaviside step function of u>, A is the reorganization energy, and u>d is 
the Debye frequency. The inverse of up represents the characteristic time scale of bath 
relaxation, and the quantum coherence can be largely preserved as ud decreases. To reduce 
the computation cost in the hierarchy equation, the high-temperature approximation is 



applied to the time correlation function, resulting in 



g(t) 



2A 

(3uj d 



-ui D \t\ 



+ i Sign(t)A- 



-U) D \t\ 



(38) 



where Sign(t) is the sign function of t. In the short-time limit (|t| — > 0), g(t) asymptotically 
follows g(t) ~ \t 2 /(3 + i\t, which is applied in the study of electron transfer 46N50|. In the 
long-time limit (|t| — > oo), the asymptotic time dependence becomes g(t) ~ 2\\t\/f3u}£> + 
i Sign(t)A/Vo, which resembles the result in the HSR model. To be consistent, the high 
temperature approximation is also used in our higher-order QKE approach. 

The parameters of our numerical calculation, e\2 = 100 cm -1 , cu^ 1 = 100 fs and T = 300 
K, are taken from Ref. 4l|. Two different site-site couplings, Jyi = 20 and 100 cm -1 , are 
chosen to test the reliability of the higher-order QKE. The results of the effective forward 
rate (/c^d) from the donor to the acceptor from the higher-order QKE and the hierarchy 
equation are plotted in Fig. |5] For the small site-site coupling (J = 20 cm -1 ), with the 
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change of the reorganization energy A, a difference up to < 40% can be resolved between 
kA<-D calculated from the hierarchy equation and from the second-order kinetic rate k^ (i.e., 
the Forster rate). As a comparison, kA^D ~ & + k^ after the leading order correction 
converges to the result of the hierarchy equation. With the Pade approximation 4g[ l^ . 
we apply the partial resummation technique, k^D ~ [k^] 2 /[k^ — k^]. This improved 
prediction agrees perfectly with the result of hierarchy equation for an arbitrary A. To 
further demonstrate that the higher-order QKE is not limited in the regime of small site-site 
coupling, we test a much larger value, J12 = 100 cm" 1 , where the transportation does not 
follow the simple hopping picture and the Forster rate k^ can be three times larger than 
the exact result. Although k^ causes an unphysical over-correction to k^ 2 \ the prediction 
after the Pade approximation compares very well with the result of the hierarchy equation in 
the whole regime of A, especially in both coherent (A < 20 cm" 1 ) and incoherent limits (A ~ 
1000 cm" 1 ). By combining the higher-order QKE method with the Pade approximation, 
the higher-order QKE can thus improve the theoretical prediction of quantum dissipative 
dynamics with a tolerable increased computation cost. 



C. Long-Range Energy Transfer in a Three-Site Bridge System 

In one model Hamiltonian of the seven-site FMO light-harvesting protein complex 
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53|, the first energy transfer pathway (sites 1 — > 2 — > 3) carries a barrier crossing 
event at site 2, and sites 1 and 3 are weakly coupled to each due to their long distance. In 
classical kinetics, such a pathway is hindered by the barrier crossing from site 1 to site 2, 
becoming less favorable compared to the alternative downhill pathway, 6 — >■ (5,7) — > 4 — > 3. 



However, our previous quantum-classical comparison [ll| has showed that the first pathway 
can dominate even at the room temperature (T = 300 K) when the electronic excitation is 
initialized at site 1. The adjustment of the energy transfer pathway is mainly caused by the 
direct energy transfer from site 1 to site 3 through multi-site quantum coherence. 

To demonstrate the long-range energy transfer phenomenon in a simple but transparent 
manner, we select the sub-system of the first energy transfer path in our seven-site FMO 
model. We further set zero dipole-dipole coupling between site 1 and site 3 (see Fig. Hb) to 
neglect the irrelevant third-order correction but focus on the leading-order corrections: 
the multi-site quantum coherence -^coherence anc ^ the bath relaxation K^ th . Following our 
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previous papers [9|, [llj, the Hamiltonian of our three-site system is given by 



H 



s 



i 280 -106 » 

-106 420 28 
i 28 j 



cm" 1 . (39) 



The Debye bath with the physiological condition is considered: A = 35 cm -1 , u^ 1 = 50 fs, 
and T = 300 K, together with the zero spatial correlation, c mn = S mn . 

The population dynamics predicted by the higher-order QKE are plotted in Fig. [6J to- 
gether with the result of the hierarchy equation. We find that in the second order, the 
quantum kinetic equation using the Forster rate is unable to reliably predict both the short- 
time quantum oscillations and the long-time kinetics. Following the Pade approximation, 
the fourth-order corrections K^ th and -^coherence are included in the rate matrix. With these 
leading-order corrections, the prediction of the higher-order QKE is significantly improved, 
compared with that of the hierarchy equation. To avoid the overcorrection of these two 
terms, the Pade approximation is used for each correction term in our calculation. Here we 
find that K^ th can quantitatively describe the profile of slow-down (< 200 fs) in the time 
evolution of populations at sites 1 and 2, although the exact quantum dynamics behaves 
as an under-damped oscillator. Further improvement requires the non-Markovian form of 
the time-nonlocal rate kernel K. Our comparison also determines that the bath relaxation 
does not affect the long-time dynamics (> 200 fs), possibly due to the fact that the bath 
relaxation time (50 fs) is still much shorter than the overall energy transfer time (~ ps). 
The more relevant multi-site coherence correction ^coherence i s shown to describe long-time 
population dynamics in a quantitatively reliable way. We find that population accumulation 
at the trap site 3 is doubled compared to the prediction using the Forster rate in 2 ps. More 
importantly, we observe a direct evidence of the long-range energy transfer: The majority 
of the fast increase in Pz(t) after t > 200 fs arises from the decrease of P\{t) instead of 
P2(t), which is consistent with the calculation of the flux network in the seven-site FMO 



model 



111 ] . Overall, the nontrivial quantum coherent effect and the bath relaxation effect 



are identified and distinguished through K^ th an< ^ -^coLr. 
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D. Quantum Phase Interference in a Closed Three-Site Loop Model 



In additional to quantum tunneling, another unique and nontrivial quantum effect is the 
interference of quantum phases. In this subsection, we use a closed three-site loop model 
as our last example (see Fig. Hb) to demonstrate the quantum phase interference successful 
predicted by the higher-order QKE. The result of the classical white noise in Eq. ( 136|) is 
extended to the quantum Debye noise, and the fourth-order corrections are also included as 
a comparison. 

For simplicity, we consider a degenerate three-site system with the following Hamiltonian, 



He 



U Jl2 ^13 
J*2 ^23 
Y "^13 ^23 o / 



(40) 



All the site-site couplings are further assumed to be the same in the amplitude, | J12I = 
I <^23 1 = I >^i3 1 = 20 cm" 1 . We assume all the other couplings are real positive and check two 
phases for the coupling between sites 1 and 3: a) J13 = 20 cm -1 , and b) J13 = 20z cm -1 . 
The imaginary phase in the second case might be generated by a coherent laser pulse. As 



shown in the study of the HSR model 29|, the quantum phase interference can cause a 
significant difference in quantum dynamics, leading to the optimization of energy transfer 
with the variation of quantum phase. 

Here we apply the Debye spectral density (A = 50 cm -1 and cu^ 1 = 10 fs) together with 
a high temperature (T = 300 K) approximation to model the bath. The initial system 
is populated at site 1. Under this particular reorganization (A = 2.5|J|), it is expected 
that the short-time quantum oscillation is suppressed. However, the nontrivial interference 
effect of quantum phase is still crucial for the incoherent dynamics. The numerical results 
of population dynamics using the higher-order QKE and the hierarchy equation are plotted 
in Fig. [7] for both conditions of J13. In the second-order, Fermi's golden rule rate cannot 
distinguish the phase of J13 and predicts the exactly same time evolution of population at 
sites 2 and 3. As shown in Fig. the higher-order QKE clarifies the effect of quantum phase 
interference. Similarly, the Pade approximation is applied to every higher-order correction 
term. For the real positive value of J13 in Fig. [Tki, is nearly negligible and the dynamics 
of the sites 2 and 3 remains degenerate, P%{t) = Pz(t), after including K^> and in the 
quantum kinetic equation. For the the imaginary value of J13 in Fig. [7b, the third-order 
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correction, /C^, causes a significant change in dynamics: 1) the population transfer out of 
the initial site 1 is accelerated; 2) the increase of Ps(t) is much faster than P2(t) in the short 
time regime; 3) a short-time uni-directional energy transfer, 1 — > 3 — > 2, is determined. 
The above phenomena can be attributed to the constructive interference for site 3 and the 
destructive interference for site 2. If site 3 is connected to a population sinker, the imaginary 
coupling of J13 can yield a hi gher energy transfer efficiency, implying the optimization on 



quantum phase accumulation 29]. The fourth-order correction in this second condition is 
much less relevant. In addition, the predictions of the higher-order QKE agrees well with 
the results of the hierarchy equation for both conditions, which again confirms the reliability 
of our methodology. 

VI. SUMMARY AND DISCUSSIONS 
A. Summary 

In this paper, we have applied a new formalism to derive a higher-order quantum kinetic 
expansion (QKE) approach to study quantum dissipative dynamics for a multi-'site' system. 
After the integration of quantum coherence and the average over the local equilibrium bath, 
we derived a closed non-Markovian quantum kinetic equation to describe the time evolution 
of system populations. In this time-convolution equation, the kinetic rate kernel K, is rigor- 
ously and systematically expanded in the order of the site-site coupling J, i.e., off-diagonal 
elements of the system Hamiltonian. The second-order rate kernel K.^ recovers the result of 
the NIBA method, and its time integration gives Fermi's golden rule rate. The higher-order 
corrections, K,^ k>2 \ include the contribution from the multi-site quantum coherence (the 
direct coherence- coherence transition) and the bath relaxation (the system-bath entangled 
population state). For a harmonic bath, the analytical expression of the kinetic kernel K 
is obtained using the displacement operator and the Gaussian cumulant expansion. Our 
higher-order QKE approach is examined in four model systems to demonstrate its reliabil- 
ity. In the Haken-Strobl-Reineker (HSR) model, the higher-order QKE leads to the identical 



kinetic mapping previously derived by the stationary approximation of coherence [29j. Un- 
der a quantum Debye noise, the prediction of the higher-order QKE together with the Pade 
approximation agrees very well with the exact result of the hierarchy equation. 
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Compared to many other theoretical approaches of quantum dissipative dynamics, the 
higher-order QKE can quantify higher-order corrections to the second-order prediction, i.e., 
the Fermi's golden-rule expansion. For example, the bath relaxation can slow down the direct 
transfer from donor to acceptor, and the exact quantum transfer rate is consistently smaller 
than the Forster rate (see Fig. For the three-site bridge model in Sec. IV C| the higher- 
order QKE predicts that the bath relaxation slows down the short-time dynamics (< 200 fs) 
whereas the multi-site coherence speeds up the long-range transfer process afterwards (see 
Fig. [6]). For the closed three-site loop model in Sec. IV Dj the quantum interference described 
by breaks the symmetry in the 'classically' incoherent dynamics (see Fig. [7]). All these 
examples have confirmed that our higher-order QKE can clarify distinct nontrivial effects 
of multi-site quantum coherence and bath relaxation, which are crucial for understanding 
nontrivial quantum effects. 



The theoretical studies in this paper and our previous two papers [111 . |29[ form a self- 
consistent methodology of quantum-classical comparison and kinetic mapping for qua ntum 
dissipative dynamics. Compared to kinetic mapping of the HSR model in Ref. 29[, the 
expansion technique in this paper has been extended from a classically white noise to an 
arbitrary quantum bath. As a result, the bath relaxation and the multi-site coherence are 
unified in the same theoretical framework. Consistent with the quantum-classical comparison 



in Ref. 



11] . the long-range energy transfer is now quantified in the detailed time evolution 



and isolated from the short-time bath relaxation. In principle, the higher-order QKE can 
be applied to an arbitrary quantum dissipative dynamic system for the investigation of the 
nontrivial quantum effects. 



B. Discussions 



The calculation of the four examples in this paper shows the numerical accuracy from the 
higher-order QKE. Using the hierarchy equation as the benchmark for the quantum Debye 
noise, the higher-order QKE can provide reliable, sometime accurate, description for both the 
average transfer rate and the detailed time evolution. The quantum dynamics of the iV-site 
system is described by the time evolution of the reduced density matrix in the N 2 -D Liouville 
space. The dimensionality of the rate matrix for the hierarchy equation grows roughly 

]\f h + 2 with the hierarchy expansion order h under the high-temperature approximation 
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for the Debye noise. On the other hand, the rate kernel K, in the higher-order QKE is always 
restricted in the N-D population subspace. Although the computation time in our method 
increases with the expansion order similarly, the Markovian approximation dramatically 
reduces the cost by changing the time-nonlocal kernels into the average rate matrix. The 
Markovian higher-order QKE can predict the overall features of time evolution. In addition, 
the partial re-summation technique based on the Pade approximation further accelerates the 
convergence of the rate matrix K^ k \ For example, with a large site-site coupling ( J12 = |^i2|), 
the re-summation of the leading-order correction, has already resulted in an almost 

quantitative prediction of the transfer rate. Thus, the higher-order QKE promises the 
potential of predicting the quantum dissipative dynamics with an acceptable computational 
cost. 

In the higher-order QKE, the quantum dynamic system is defined as a general iV-'site' 
network form. The so-called 'site' can be further generalized as any basis set of the system 
Hilbert space. Thus, the higher-order QKE is not restricted in energy transfer and electron 
transfer, but can extended to other quantum dynamic processes. The surrounding bath is 
also defined generally in the higher-order QKE. Many sophisticated deterministic or stochas- 
tic methods, such as the hierarchy equation, the SC-IVR, the QUAPI, the polaron-based 
methods, the path integral Monte Carlo, etc., are based on a presumed harmonic bath, which 
is in general a good approximation under many conditions. However, theoretical methods 
for the anharmonic environments is also highly required. Since the expressions of the time- 
nonlocal rate kernels KP^ in Eqs. ( |T4l) and ( TT5l) are independent of the bath model, the 
formulation in the Hilbert space, e.g, Eq. ( 12^1) . can work as the starting point of studying 
quantum dissipative dynamics in a complicated environment. Numerical implementation 
will be worth in the higher-order QKE along this direction. 

Our current derivation needs an incoherent preparation for the initial product state, 
which can be a strong assumption considering the experimental designs, e.g., a coherent laser 
pulse can generate different initial states. With an additional expansion, the improvement 
of including the initial quantum coherence can be derived in the future. Together, the time 
evolution of quantum coherence needs to be derived after the higher-order QKE. A more 
important question is about the expansion parameter of the higher-order QKE, i.e., the 
site-site coupling J mn . The different site-site couplings are also not necessarily on the same 
order of magnitude. To solve this difficulty, we can introduce the sub-system concept and 
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construct the expansion based on the weak coupling between sub-systems. For exam ple, 
the multichromophore Forster resonance energy transfer (MC-FRET) rate theory 54m 6 1 
can serve as the second-order prediction in the extended higher-order QKE of multiple sub- 
systems instead of multiple sites. The higher-order corrections then help reveal nontrivial 
quantum effects beyond the MC-FRET result. Overall, the higher-order QKE requires future 
improvements to solidify its construction and extend its applications. 
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FIG. 1: The diagrammatic representation of dynamic trajectories in the second- (a) and third- 
order (b) quantum rate kernels. Each circle with a single letter denotes a decoupled population 
state, e.g, P n (t)pj^; each circle with two letters denotes a system-bath entangled coherence state, 
e.g., pmnif) = Ylbv Pmb,nb' {b' \. The duration time, e.g., ri^,..., spanned at each density state 
(population or coherence) is given beneath its circle. Each dashed line between a pair of single- 
letter circles represents a nonzero interaction J between these two sites in H$- Each directed 
curve represents a transition from one density state to the subsequent one, where the inducing 
interaction, e.g., J mn , is provided. Each integrated diagram composed of circles and connected 
curves describes one trajectory in the quantum rate kernels: a) a typical term in K,^ of Eq. (|29p 
and b) a typical term in K,^ of Eq. (|30|) , 
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a) b) 




Ti X 3 



FIG. 2: Three typical dynamic trajectories of density states in the fourth-order quantum rate kernel, 
^bath-mn" Each diagram represents a distinct behavior of the bath relaxation in Eq. f)31 1) (see text). 
Here all the symbols, circles and lines have been explained explicitly in Fig.[H except for the dashed 
circles that represent system-bath entangled population states, e.g., p nn (t) = Ylbb' Pnb,nb'(t)\b)(b'\. 
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FIG. 3: Four typical dynamic trajectories of density states in the fourth-order quantum rate kernel, 
^coherence-ran" Each diagram represents a distinct behavior of the multi-site quantum coherence in 
Eq. (|32p (see text). Here all the symbols, circles and lines have been explained explicitly in Fig. [U 
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a) 




FIG. 4: The schematic Hamiltonian diagrams of the three example quantum networks studied in 
Sec.|V] a) the two-site system, b) the bridged three-site system, and c) the closed-looped three-site 
system. Here each circle with a number represents a 'site' (basis) of the system. Each dashed line 
denotes a nonzero site-site coupling J mn . The height of each circle denotes the relative energy e m at 
each site. The closed-looped three-site system in c) actually forms a triangle network geometrically. 
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FIG. 5: The effective forward rate (kAi-D) from the donor to the acceptor in the two-site sys- 
tem (Fig. Hk) calculated using the higher-order QKE and the hierarchy equation. The detailed 
parameters are provided in Sec. IV Bl Among them, the site-site coupling is different in the two 
panels: J = 20 cm -1 in a) and J = 100 cur 1 in b). The dashed line denotes the Forster rate (i.e., 
the second-order rate of QKE). Both the dot-dashed and the solid lines include the fourth-order 
correction of bath relaxation, whereas the solid lines include the additional Pade approximation 
(see text). As a comparison, the data from the hierarchy equation are plotted as the solid dots. 
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t(ps) 

FIG. 6: The population dynamics for the three-site bridge model (Fig. 0b) with the Hamiltonian 
in Eq. f|39|) calculated by the higher-order QKE with three different rate matrices and by the 
hierarchy equation, respectively. The bath parameters are provided in Sec. IV CI From the top to 
the bottom, three distinct sets of curves represent the time evolution of Pi(t), P2(t), and Ps(t). 
Here the dotted lines are the results from K^; the dashed lines are the results from + K^.,; 
the solid lines are results from + -f^ath + -^coherence- ^ ne P a de approximation is applied to all 
the fourth-order corrections. As a comparison, the results of the hierarchy equation are plotted in 
the solid lines highlighted by the solid dots. 
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FIG. 7: The population dynamics for the closed-looped three-site model (Fig. Ht) with the Hamil- 
tonian in Eq. (I40p calculated by the higher-order QKE with rate matrices of three orders and by 
the hierarchy equation. The bath parameters are provided in Sec. IV Dl The left panel a) presents 
the results of J13 = 20 cm -1 ; the right panel b) presents the results of J13 = 20i cm" 1 . In each 
panel, the time evolution curves are labeled by -Pi, 2, 3 f° r the three sites. For the higher-order 
QKE, the dotted lines are the results from K^; the dashed lines are the results from + K^; 
the solid lines are results from + + . The Pade approximation is applied to all the 
higher-order corrections. As a comparison, the results of the hierarchy equation are highlighted by 
symbols (circles, diamonds, and rectangles for Pi, P2, and P3, respectively). 
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